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This paper puts forth two new closure models for the proper orthogonal decomposi- 
tion reduced-order modeling of structurally dominated turbulent flows: the dynamic 
subgrid-scale model and the variational multiscale model. These models, which are con- 
sidered state-of-the-art in large eddy simulation, together with the mixing length and the 
Smagorinsky closure models, are tested in the numerical simulation of a 3D turbulent 
flow around a circular cylinder at Re = 1, 000. Two criteria are used in judging the perfor- 
mance of the proper orthogonal decomposition reduced-order models: the kinetic energy 
spectrum and the time evolution of the POD coefficients. All the numerical results are 
benchmarked against a direct numerical simulation. Based on these numerical results, we 
conclude that the dynamic subgrid-scale and the variational multiscale models perform 
best. 
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1. Introduction 

Reduced- order models (ROMs) of structurally dominated turbulent flows are central to 
many applications in science and engineering, such as fluid flow control (see for example 

19991 ICohen et aZ.||2003 



Ito & Ravindran 1998 Graham et al 



Lehmann et aL||2005[ |Hoepffner et aZ.||2006||Bagheri et al.|2009 



Bcrgmann et aZ.||2005 



Barbagallo et al.j 2009 



Ahuja &: Rowley||2010| |Akhtar fc Nayfeh||2010[) and da ta assimilation of atmospheric 



and oceanic flows (Luo et al. 2007 Daescu & Navon 2008 Fang et al. 2009). Both 



computational efficiency and physical accuracy are needed for the success of these ROMs 
in practical applications. Striking a balance between efficiency and accuracy in ROMs 
of turbulent flows is, of course, challenging. Indeed, it is clear that the fewer the modes 
retained in the ROM, the more efficient the ROM is. Preserving the physical accuracy of 
the resulting ROM, however, becomes challenging, since the modes that are not retained 
in the ROM representation of the underlying turbulent flow need to be modeled. The 
quest of balancing the computational efficiency and physical accuracy represents one of 
the main challenges in ROMs for turbulent flows. 

One of the most successful ROM strategies for structurally dominated turbulent flows 
has been the Proper Orthogonal Decomposition (POD) (see for example Holmes et al. 



1996 Sirovich 1987). POD starts with data from an accurate numerical simulation (or 
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physical experiment), extracts the most energetic modes in the system, and utilizes a 
Galerkin procedure that yields a ROM of the underlying turbulent flow. The first proper 
orthogonal decomposition reduced-order model (POD-ROM) for the turbulent boundary 
layer was proposed in Aubry et al. (1988). This model truncated the POD basis and 



used an eddy viscosity-based approximation to model the effect of the discarded POD 
modes on the POD modes kept in the model. This POD-ROM yielded good qualitative 
results, considering the coarseness of the approximation. The criterion used to assess 
the accuracy of the model was the intcrmittency of bursting events in the turbulent 
boundary layer. This POD-ROM was further investigated numerically in two subsequent 



papers (Podvin & Lumley 1998 



Podvin 2001). The model reproduced the qualitative 

by adding new POD modes 



physics of the turbulent boundary layer well. Furthermore, 
to the model, the accuracy of the model was increased. 

Despite their initial success, POD-ROMs have generally been limited to laminar flows 
and relatively few reports on closure modeling strategies for turbulent flows have ap- 



peared in the literature i 


Aubry et al. 1988 


Podvin & Lumley 1998 


Podvin 2001 


Rcmpfer 


fe Fasel||1994[ |Rempfer| 
Karniadakis||2004| |Buffc 


1996;|Cazemier et al. 1998|IMa & Karniadakis||2002[ |Sirisup & 
mi et a/.||2006| INoack et al.\ 20021 120031 120051 120081 lUllmann & 



in traditional turbulence modeling, such as large eddy simulation (LES), where liter- 
ally hundreds of closure models have been proposed and investigated (see for example 
Sagaut 2006) over the same time period. This disparity in closure modeling between 



POD reduced-order modeling and classical turbulence modeling seems even more dra- 
matic considering that the concept of an energy cascade, which is a fundamental modeling 
principle in LES, is also valid in a POD setting. Indeed, the validity of the extension of 
the energy cascade concept to the POD setting was investigated numerically in [Couplet] 
et al. (2003). The authors have investigated the energy transfer among POD modes in a 



non-homogeneous computational setting. By monitoring the triad interactions due to the 
nonlinear term in the Navier-Stokes equations, they have concluded that the transfer of 
energy among the POD modes is similar to the transfer of energy among Fourier modes. 
Specifically, they found that there is a net forward energy transfer from low index POD 
modes to higher index POD modes and that this transfer of energy is local in nature 
(that is, energy is mainly transferred among POD modes whose indices are close to one 
another) . This study (see also Noack et al. 2002 ) clearly suggests that LES ideas based 
on the energy cascade concept could also be used in devising POD-ROMs. 

One of the main reasons for the scarcity of closure models for POD-ROMs of turbulent 
flows is the impractical cost of standard LES closure models employed in a POD-ROM 
setting. Indeed, most of the computational cost of a POD-ROM lies in assembling the 
vectors, matrices and tensors of the ROM. This, however, is hardly a problem for POD- 
ROM, since the vectors, matrices and tensors are assembled only once, at the beginning of 
the POD- ROM simulation, and reused at every time step. Standard (nonlinear) LES clo- 
sure models, however, introduce new vectors and matrices that need to be recomputed at 
every time step. Thus, a straightforward numerical discretization of such closure models 
would come at a huge computational cost, rendering the resulting POD-ROMs imprac- 
tical. 

In the past few years, a number of strategies have been introduced to treat nonlinear 
terms in POD-ROMs. These include interpolatory methods such as the empirical inter- 
polation method ( jBarrault e^"Z1|2004| |Chaturantabut et a/.||2010| |Galbally et aL||2010[ ), 
the closely related group finite element approach ( Dickinson fc Singler|2010 ) and a novel 
two- level discretization method ( Wang et aL|201l" ). The latter approach is best suited for 
this study since it does not constrain the nonlinear term to lie within a predefined set. 
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This approach is based on a two-level discretization of the vectors, matrices and tenstors 
of the POD-ROM, in which all the terms are computed on the fine grid, except for the 
nonlinear closure model terms, which are computed on a coarser grid. In Wang et al. 
(2011), numerical simulations of a turbulent flow past a 3D cylinder at Re = 1,000 with 
a standard LES closure model (Smagorinsky 19631 have shown that the new two-level 
discretization is both computationally efficient and physically accurate. Indeed, the new 
two-level algorithm decreased by more than an order of magnitude the CPU time of the 
standard one-level algorithm, without compromising the physical accuracy. 

In this report, we use the two-level algorithm proposed in Wang et al. (2011) to dis- 
cretize two new POD-ROMs, inspired from state-of-the-art LES closure modeling strate- 



gies: the dynamic subgrid- scale (DS) model (Germano et al. 1991 Meneveau et «Z.||1996 



Porte- Agel et «L|2000 ) and the variational multiscale (VMS) model ( Hughes et aL|2000 ). 



We also consider the standard mixing-length closure model proposed in Aubry et al. 



2002 



(1988) and the Smagorinsky model proposed in [Wang et aZ.| (2011 ) (see also Noack et al. 



Ullmann & Lang 2010), both being standard LES closure models. All four POD- 



ROMs are tested in the numerical simulation of a 3D turbulent flow around a circular 
cylinder at Re — 1, 000. Two criteria are used in judging the performance of the POD- 
ROMs: the kinetic energy spectrum and the time evolution of the POD coefficients. All 
the numerical results are benchmarked against a direct numerical simulation. 

The rest of the paper is organized as follows: The general methodology used in the de- 
velopment of POD-ROMs is presented in §[2j The four POD closure models are described 
in §[3] and are investigated numerically in §[4] Finally, conclusions and several research 
directions currently pursued by our group are provided in §[5j 



2. POD Reduced-Order Modeling 

We now present the general approach used in the development of POD-ROMs. We 
start by briefly describing the POD methodology. For more details, the reader is referred 
to Sirovich (19871; Holmes et al. (1996). To this end, we consider the numerical solution 



of the incompressible Navier-Stokes equations (NSE): 

u t - Re -1 Au + (u ■ V)u + Vp - 

V • u = 



0, 



(2.1) 



where u is the velocity, p the pressure and Re the Reynolds number. The POD basis 



is generated by post-processing typical data from the numerical simulation of (2.1). If 
y = {y(-,i) SK | t G (0, T)} (with T~L a Hilbert space) represents a simulation of the 
NSE, then the first POD basis vector is the function that maximizes the time-averaged 
projection of y onto itself, 



1 

max — 

<p^:'H,\\ip\\-H = l 1 



I <y(-> *),*>(•)) 



Hi 



dt. 



(2.2) 



Subsequent vectors, (p k , are determined by seeking the above maximum in the orthogonal 
complement to 



X fc - X =span{ ¥ > 1 ,...,v> fe _ 1 }, 2^k^N, 



H, 



(2.3) 



where N is the rank of y. If we choose T-L = £2 and y represents a single simulation, the 
POD basis functions satisfy the Fredholm integral equation 



/ R(x,xV t (xVx' = A^(x), 



(2.4) 
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where 



R(x,0 = ~J y(x,f)y*(x', 



t) dt 



(2.5) 



is the spatial autocorrelation kernel. There are natural extensions of this definition that 
accommodate multiple simulations. In practice, either the time average of each simulation 
or the steady state solution is removed, so that y contains fluctuation from the mean 
(or a centering trajectory), e.g., y(x, f) = u(x, t) — U(x) ( Holmes et aZ.|1996 1. Note that 
each POD basis vector ip k represents a weighted time average of the data y. Thus, these 
basis vectors preserve linear properties (such as the divergence- free property). 

A POD basis enables a reduced representation of the simulated data, and thus can 
be viewed as a compression algorithm. Utilizing the POD basis to obtain efficient ap- 



proximations to (2.1 1 is achieved using the POD basis in a Galerkin approximation, and 



employing the fact that the POD basis vectors are mutually orthogonal. A POD-ROM 
of the flow is constructed from the POD basis by writing 



u(x, t) rj u r (x, t) 



U(x)+^ aj W Vj (x), 



(2.6) 



where U(x) is the centering trajectory, {<£j} r j=i are the first r POD basis vectors, and 
{aj{t)Yj=i are the sought time- varying coefficients that represent the POD-Galerkin tra- 
jectories. We now replace the velocity u with u r in the NSE (2.1 1, and then project 
the resulting equations onto the subspace X r . Using the boundary conditions and the 
fact that all modes are solenoidal, one obtains the POD Galerkin reduced-order model 
(POD-G-ROM): 



~dt 



ip + ((u r • V)u r ,<p) 



Re 



D(u r ),Vv? =0 V93 e x r 



(2.7) 



where D(u r ) := (Vu r + (Vu r ) )/2 is the deformation tensor of u r . We note that, since 



the computational domain that we consider is large enough, the pressure terms in (2.7 1 



can be neglected (for details, see |Noack et aZ.||2005| |Akhtar et al.||2009[ ). The POD-G- 
ROM (2.7| yields the following autonomous dynamical system for the vector of time 



coefficients, a(£): 



b + Aa + a T Ba, 



(2.8) 



where b, A, and B correspond to the constant, linear, and quadratic terms in the nu- 



merical discretization of the NSE (2.1 1, respectively. The initial conditions are obtained 
by projection: 



o i (0) = < Vi ,u(. J 0)-U(-))«, j = l,. 



. ,r. 



(2.9) 



The finite dimensional system (2.8) can be written componentwise as follows: For all 
fc = l,...,r, 



a>k(t) = h + 2J 4m« m (i) + X! X! B km n a n (t)a m {t), 

in— 1 m— 1 n— 1 



(2.10) 
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^^-(^U.VU)-^^^^), (2 . n) 

A km = -(<p h ,U- VcpJ - (<p k ,<p m • VU) - A fv<p k , V ^» + V( ^ , (2 . 12) 

B fcm „ = -(Vfe, Vm • Vv„). (2.13) 



3. POD Closure Models 

In this section, we present the four POD closure models investigated numerically in §[4j 
To this end, we start by describing the filtering operation utilized and the spatial length- 
scale S used in the POD closure models. Both are needed in order to define meaningful 
LES-inspired POD closure models. 

3.1. POD Filter 

In LES, the filter is the central tool used to obtain simplified mathematical models that 
are computationally tractable. The filtering operation is effected by convolution of flow 
variables with a rapidly decaying spatial filter g$, where S is the radius of the spatial filter. 
In POD, however, there is no explicit spatial filter used. Thus, in order to develop LES- 
type POD closure models, a POD filter needs to be introduced. Given the hierarchical 
nature of the POD basis, a natural such filter appears to be the Galerkin projection. For 
all u £ X, the Galerkin projection u g X r is the solution of the following equation: 

(u-u,v>) = V¥>eX r . (3.1) 



The Galerkin projection defined in (3.1 ) will be the filter used in all POD closure models 
studied in this report. 

3.2. POD Lengths cale 

Next, we introduce the lengthscale 5 used in the POD closure models. We emphasize that 
this choice is one of the fundamental issues in making a connection with LES. Indeed, 
we need such a lengthscale (5) in order to define dimensionally sound POD models of 
LES flavor. 



To derive the lengthscale <5, we use dimensional analysis. Aubry et al. (1988) defined 
/>, a dimensionally sound lengthscale for a turbulent pipe flow. In fact, this lengthscale 
was only defined implicitly, through the turbulent eddy viscosity vt := it> Z>. Indeed, 



equation (22) in Aubry et al. ( 1988 ) reads 



v T := U> Z> = J ° W> l> > \ 1/2 , (3.2) 



' X 2 J 2 («i> j u i>d ) dx 2 



where repeated indices denote summation, the subscript > denotes unresolved POD 
modes, 



(/) = — L- / / f(x,t)d Xl dx 3 (3.3) 



denotes the spatial average of / in the homogeneous directions (here x\ and X3), and 
L\, Li and X 2 are the streamwise, spanwise, and wall- normal dimensions of the computa- 
tional domain, respectively. Note that the authors only considered the wall region, not the 
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N 

entire pipe flow. In (3.2), the following notation was used: = a* ipj, Uiy Ui> 

j=r+l 

3 



y"\i> Mj>, and u i>t j 



dui 



dx-j 



Note that a quick dimensional analysis shows that the 



i=l 

quantity defined in (3.2) has the dimensions of a viscosity. Indeed, 

,2 



[u T ] 



S S 



1/2 



m 
s 



(3.4) 



In Appendix B of Aubry et al.\ ( |1988[ ), the authors have further simplified (3.2) and 
expressed vt in terms of the first neglected POD modes: 



Z^(k,n) A k 



Xi Li L3 



(k,ra) \ 



(n) 



f,f 2 £)$J n W"^ dx 2 -k\-k 

JO »k »k ^ 1 



'-.(") rW n ) 



where the triplets (k, n) are the first neglected POD modes. 



In equation (9.90) in Holmes et al. (1996), the authors define another dimensionally 



sound turbulent viscosity 



vt '■= it> l> = 



X-2 



X 2 Jo (u i>d u 4>J ) 1/2 



dxo 



(3.6) 



A quick dimensional analysis shows that the quantity defined in (3.6) also has the di- 
mensions of a viscosity. 

We can use the two definitions of vt in (3.2) and (3.6) to define a lengthscale l>. We 
obtain 



and 



respectively. 



J* 2 (ui> Ui>) dx 2 
X 2 f 2 (Ui>,j u i>t j) dx 2 



1/2 



X 2 J() (UiyjUiyj) 



(3.7) 



(3.8) 



For our 3D flow past a cylinder, both (3.7) and (3.8) are valid candidates for the 
definition of the lengthscale S. The only modification we need to make (due to our 
computational domain) is to replace the horizontal averaging by spanwise averaging and 
take double integrals in the remaining directions. Specifically, we have 



6 := 



and 



Iq 1 Io 2 ( u i> u i>)dxidx 2 \ / 
< Jo' Jo 2 ( Ui > <J Ul > .3) dxi dx2 J 

" L2 (Ui>Ui>) 



L\ L 2 j 



1/2 



( u i>,j u i>,j) 



dxi dx 2 



(3.9) 



(3.10) 
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3.3. POD Closure Models 



We are now ready to present the four POD closure models that will be investigated 
numerically in §[4j 



The POD-G-ROM (2.7 1 can be used for laminar flows. For structurally dominated 
turbulent flows, however, the POD-G-ROM simply fails ( Wang et aL|20li" ). The reason 
is that the effect of the discarded POD modes {tp r+1 , . . . , (p N } needs to be included 
in the model. For turbulent flows, the most natural way to tackle this POD closure 
problem is by using the eddy viscosity (EV) concept. Indeed, most closure models used 
in turbulence modeling are based on this EV concept, which states that the role of the 
discarded modes is to extract energy from the system. The concept of energy cascade, 
which is well established in a Fourier setting, has been recently confirmed in a POD 



setting in the numerical investigations in Couplet et al. (2003). Thus, using LES inspired 



EV closure models in POD-ROM represents a natural step. 

In this section, we propose two new POD closure models: the dynamic subgrid-scale 
model and the variational multiscale model. We emphasize that, although these models 



were announced in Borggaard et al. ( 2008 ) , this study represents their first careful deriva- 



tion and thorough numerical investigation. We also numerically test the mixing-length 
flAubiy et aL||1988[ ) and Smagorinsky ( |Noack et al]\2002\ |Ullmann fc Lang||2010| |Wang| 
et aZ.||2011[ ) POD closure models. 

Since all four POD closure models are of EV type, we first present a general EV POD- 
ROM framework. Then, for each closure model, we specify the changes that need to be 
made to this general framework. The general EV POD-ROM framework can be written 
as: 



a = b + b(a) + A + A(a) a 



a T Ba, 



(3.11) 



which is just a slight modification of the POD-G-ROM (2.8). The new terms in (3.11) 



(the vector b(a) and the matrix A (a)) correspond to the numerical discretization of the 
POD closure model. In componentwise form, equation (3.11) can be written as 



flfc(i) = (b k + b k (a)j + ^2 { A km + A km (a)j a m (t) 

m—1 

r r 

+ ^ ^2 B k mn a n {t)a m {t), (3-12) 



where b k , A km , and B kmn are the same as those in equations (2.8) and 6fc(a) and Afc m (a) 
depend on the specific closure model used. 



3.3.1. The mixing-length POD reduced-order model (ML-POD-ROM) 

The first POD closure model was the mixing-length model proposed in |Aubry et al.\ 
(1988). This closure model is of EV type and amounts to increasing the viscosity coeffi- 
cient v by 

vml = av T = a U ML L ML , (3.13) 



where Uml and Lml are characteristic velocity and length scales for the unresolved 
scales, and a is an 0(1) nondimensional parameter that characterizes the energy being 
dissipated. Using the EV ansatz in (3.13), the mixing-length POD reduced-order model 
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(ML-POD-ROM) has the form in fl3.1ip, where 



&fc(a) = -"ML V<p 



vu + VU T 



-Afcm(a) = —VML 



(3.14) 
(3.15) 



The parameter a is expected to vary in a real turbulent flow, and different values of 



a may result in different dynamics of the flow (Aubry et al. 1988 Holmes et al 



et al. 



Podvin & Lumley 1998 Podvin 2001 ). There are also different ways to define vt hi ( 3.11 ) 
relatio n (|3.2 ) was used in Aubry et al. (1988), whereas relation (3.6) was used in Holmes 



1996 



(1996). We also mention that several other authors have used the ML-POD-ROM 



(3.13) (see for example Borggaard et al. 2008 Wang et al. 2011). Improvements to the 



mixing- length model (3.13) in which the EV coefficient is mode dependent were proposed 



in Rcmpfer & Fasel (1994); Cazemier et al. (1998); Podvin (2009). 



3.3.2. The Smagorinsky POD reduced-order model (S-POD-ROM) 

A potential improvement over the simplistic mixing-length hypothesis is to replace the 



constant vml in (3. 14)- (3. 15) (which is computed only once, at the beginning of the 
simulation) with a variable turbulent viscosity (which is recomputed at every time step) , 
such as that proposed in Smagorinsky ( 1963[ ). This yields a POD closure model in which 
the viscosity coefficient is increased by 



2(C s S) 2 ||D(u r 



(3.16) 



where C$ is the Smagorinsky constant, 5 is the lengthscale defined in §3.2 and ||JD(u. 



is the Frobenius norm of the deformation tensor D(u r ). Using the EV ansatz in (3.16), 



the Smagorinsky POD reduced-order model (S-POD-ROM) has the form (3.11), where 

VU T 



6 fc (a) = -2(C s Sy 
A km (a) = -2{C S S) 2 



V<p fc ,||B(u r ) 



V Vfe , ||D(u r 



VU 



(3.17) 
(3.18) 



The S-POD-ROM (3.17)-(3.18) was proposed in Borggaard et al. (2008) (see also Noack 



et al. 2002) and was used in the reduced-order modeling of structurally dominated 3D 



turbulent flows in |Wang et al. (2011); Ullmann & Lang (2010). Its advantage over the 



ML-POD-ROM ( 3.14 )-( 3.15) is obvious: the latter utilizes a constant EV coefficient at 



every time step, whereas the former recomputes the EV coefficient (which depends on 
||]D(u r )||) at every time step. To address the significant computational burden posed by 
the recalculation of the Smagorinsky EV coefficient at every time step, a novel two-level 
discretization algorithm proposed in Wang et al. (20111 is employed in §|4j 



3.3.3. The Variational Multiscale POD reduced-order model (VMS-POD-ROM) 

The VMS method, a state-of-the-art LES closure modeling strategy, was introduced 



Hughes et al. (2000 2001q|& ). The VMS method is based on the principle of locality 



of energy transfer, i.e., it uses the ansatz that energy is transfered mainly between the 



neighboring scales. In Couplet et al. (2003), the transfer of energy among POD modes 



for turbulent flow past a backward-facing step (a non-homogeneous separated flow) was 
investigated numerically. In their report, it was shown that the Fourier-decomposition 
based concepts of energy cascade and locality of energy transfer are also valid in the POD 
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context (Figures 3 and 4 in Couplet et aZ.||2003 ). Thus, VMS closure models represent a 
natural choice for POD-ROM. 

To develop the VMS POD closure model, we start by decomposing the finite set of 
POD modes X r into the direct sum of large resolved POD modes X£ and small resolved 
POD modes X r Q : 



X r = X r L 8 X r s , where 

X r L := span{ip l7 tp 2 , . . . ,<p rL } and 

X S : = S P an {V-r i +l ) ¥'r I ,+2 J ■ ■ ■ ,<Pr} • 



(3.19) 
(3.20) 
(3.21) 



Accordingly, we decompose u r into two components: representing the large resolved 
scales, and representing the small resolved scales: 



uf 1 + u£, 



(3.22) 



where 



u r = U + a 3 <p k , 

3=1 

r 



(3.23) 
(3.24) 



The two components and represent the projections of u r onto the two spaces X r L 
and X r s , respectively. The general POD- ROM framework (3.11) can now be separated 
into two equations - one for a L (the vector of POD coefficients of u^) and one for a s 
(the vector of POD coefficients of uf). The Variational Multiscale POD reduced-order 
model (VMS-POD-ROM) applies an eddy viscosity term to the small resolved scales only, 
following the principle of locality of energy transfer. The VMS-POD-ROM reads: 



a s 




h s 


+ A r 


' a L " 

a s 


+ 




+ 


' a L ' 

a s 


T r 

B 


a L " 









A s {a s ) J 



(3.25) 



The finite dimensional system (3.25) can be written componentwise as follows: 



ifc (0 = b k + J2 A krna m {t) + A % a i (*) + EE B kmn a n (t)a m (t), (3.26) 

3=1 



rn — 1 



m— 1 n— 1 



Vfc = 1 



l(t) = b s k + J2 A r km a m (t) + J2 ( A t, ■ < l )".' h 



m—l 

r r 



^ ^ Bkmn a n(t) a m(t) 



(3.27) 



Wk = r L + 1, .. . ,r, 
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where 



A r = 

km 



A L 



A s 



Akj ( a ) 



2 / VU + VU T 

-(^u.vu)--(v^ 

-( Vfc ,U-V Vm )-(v fc)¥ » m -VU), 



Re V ' rk, ~ ~T 
■(<Pk,<Pm- V Vn)> 



Re 



-2(C S <5) 2 ( Vy fc ,l|P(uf + U)|| V¥?J ' + Vy/ 



(3.28) 
(3.29) 

(3.30) 

(3.31) 
(3.32) 

(3.33) 
(3.34) 



We emphasize that the system of equations (3.25) is coupled through two terms: (i) 



a Ba, which represents the nonlincarity (u r • V)u r ; and (ii) A r a, which represents the 
term (u r • V)u r linearized around the centering trajectory U. The difference between the 
VMS-POD-ROM ( |3~25| )- p34| and the S-POD-ROM fjTl^ - flOil ) is that the former acts 
only on the small resolved scales (since the Smagorinsky EV term (Cg 8) 2 ||B(u° +U)|| is 
included only in the equation corresponding to a s ), whereas the latter acts on all (both 
large and small) resolved scales. 

The VMS-POD-ROM (|3.25|)-(|3.34|) was announced inlBorggaard et al.\ (I2008I). This 



study, however, represents its first careful derivation and through investigation in the 
numerical simulation of a 3D turbulent flow. We note that a fundamentally different 
VMS LES closure model that utilizes the NSE residual was proposed in|Bazilevs et al.\ 



(2007); this model was used in a POD setting in Bergmann et al. (20091). Yet another 



VMS-POD-ROM, inspired from the numerical stabilization methods developed in |Layton| 
pOOl; iGuermondl (fT^99b; iJohn & Kayal pOOSl; iJohn & Kindll pOlob, was proposed, 



analy z ed an d tested in |Iliescu fc Wang| ( |2010[ ). We em phasize that the VMS-P OD-ROM 
( 3.25 )-([3.34 1 is different from both the model used in Bergmann et al. (2009) and that 
used in |Iliescu fc Wang (2010). 



3.3.4. Dynamic Subgrid-Scale POD reduced-order model (DS-POD-ROM) 

For all thr ee POD-ROM c losure mode ls defined up to this poi nt (i.e ., ML- POD-ROM 
( |37l4l )- pTT5] >, S-POD-ROM pl^ - pls) ), and VMS-POD-ROM p^-plil)), the defi- 
nition has been entirely phenomenological. Indeed, arguing that the role of the discarded 
POD modes is to extract energy from the system, we used an EV ansatz to derive 
closure models of increasing complexity and physical accuracy. The dynamic subgrid- 
scale (DS) POD-ROM closure model is also of EV type. Its derivation, however, needs 
a precise definition of the filtering operation. The DS closure model has its origins in 
LES, where it is considered state-of-the-art (see for example Sagaut 2006). In LES, the 
filtering operation is effected by convolving the flow variables with a rapidly decaying 
spatial filter. In POD, the filtering operation is effected by using the POD Galerkin pro- 
jection described in ^3.1 (see (3.1)). To derive the precise POD filtered equations, we 
start with the NSE (12. 1\ in which the velocity u is replaced by its POD approximation 
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u(x, t) « u r (x,i) = U(x) + J2 V j=i a j(^)Vj( x ) m (2-6), and obtain 
du r 



dt 



- Re 1 Au r + (u r • V) u r + Vp = 0. 



(3.35) 



Using the fact that V • u r = in (3.35), we get (u r ■ V) u r = V • (u r u r ). Thus, ( |3.35 ) 
can be rewritten as 



du r 

~dt 



- Re 1 Au r + V • (u r Ur) + Vp = 0. 



(3.36) 



Applying the POD filtering operation (3.1) to (3.36), using the fact that the POD 



Galerkin projection is a linear operator, and assuming that differentiation and POD 
filtering commute, we obtain 

(yU 

-^-i? e - 1 Au r + V-(u7u r ") + Vp = 0. (3.37) 

We note that, if filtering and differentiation do not commute, one has to estimate the 
commutation error (see for example |Vasilyev fc Goldstein 2004 Vasilyev et al. 1998 



|Berselli et aL|2006 ). We also note that, since the POD filtering operation is the Galerkin 
projection (3.1 ), u r = u r . For consistency with the nonlinear term notation, we still use 



the u,. notation in what follows. 



The POD filtered equation (3.37) can be rewritten as 
du r 

Ik 

where 



i?e _1 Au r + V • (u r u r ) + V • (r. r ) + Vp = 0, 



U r U r 



U r U r 



(3.38) 



(3.39) 



is the POD subfilter-scale stress tensor. Thus, the POD-G-ROM (2.7) amounts to setting 



T r = 0. For turbulent flows, as we have already mentioned, this approximation is flawed. 
Thus, one needs to address the POD closure problem, i.e., to model the POD sufiltcr- 
scale stress tensor r r in terms of POD filtered velocity u r . We note that the POD closure 
problem is exactly the LES closure problem, in which the spatial filtering is replaced by 
POD Galerkin projection. For all three POD-ROM closure models defined so far in this 
section (i.e., ML-POD-ROM (|3T4|-((3T5l), S-POD-ROM (|3T7]-(|3T8|, and VMS-POD- 



ROM (3.25 )-(3.34 )), the closure problem has been addressed by assuming an EV ansatz 
for T r . The DS-POD-ROM employs an EV ansatz as well; specifically, the Smagorinsky 
model is utilized: 



:= (C s <5) 2 ||D(ur 



(3.40) 



in which Cs is not a constant (as in the Smagorinsky model), but a function of space 
and time, i.e., Cs = Cs(x,t). To compute Cs(x,<), we follow the LES derivation in 



Sagaut (2006) and replace the LES spatial filtering with the POD Galerkin projection. 



Since there are two spatial filters in the LES derivation of the DS model, we define a 



second POD Galerkin projection (in addition to that defined in (3.1 )): For all u G X, the 
second (test) Galerkin projection u g X fl (where R < r) is the solution of the following 
equation: 

(u-u,p)=0 V^eX 8 



Applying the second POD filtering operation (3.41) to (3.37), we obtain 



<9u,. 



Re x An r + V ■ (u r u r ) + V ■ (T r ) + Vp = 0, 



(3.41) 



(3.42) 
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where 

TV = u r u r — u r u r (3.43) 

is the second POD suhTtcr-scale stress tensor. We note that the following identity (called 
the "Germano identity" in LES) holds: 

L r - r' r . (3.44) 



u r u r — u r u, 



u r u r u r u, 



where L r = u r u r — u r u r and r r = u r u r — u r u r . We assume the same EV ansatz for 
the two POD subfhter-scale stress tensors, r r and T r : 

T r w -2(C S S) 2 ||ro(n,.)||D(n r ) (3.45) 
r r « ~2(C S S) 2 ||D(u r )||D(u r ), (3.46) 



where 5 is the filter radius used in the second POD filtering operation (3.41 ). Assuming 

(3.47) 



that Cs remains constant under the second POD filtering (3.41), we obtain: 

TV w -2 (Cs <5) 2 ljD(u r )|| D(ur) « -2 (C s S) 2 ||D(u^jlD(u r ). 
Plugging ( |3.45| and ( |3.47[ ) into ( |3.44[ ) we obtain: 



2(Cs<5n|B(u r )||B(u r ) = u r u r -u r u r -2(Cs5) 2 ||D(H r )||B(u r ). (3.48) 



We note that Cs is the only unknown in (3.48), all the other terms being computable 



quantities. Since all the terms in (3.48) are tensors, the unknown Cs cannot satisfy all 



nine equations. Thus, the following least squares approach is considered instead: 

min r , _ „ \ - — • ~ 

TV&r - u r u r ) - 2 (C s S) 2 ||B(u r )|| B(u r ) + 2 (C s S) 2 ||B(u r )|| B(n r ) 

2 (C s S) 2 ||B(u^lD(u r ) + 2 (C s S) 2 ||B(n r )|| B(n r )l (3.49) 



r 2 



The solution C s (x,t) of (3.49) is: 



Cf(x,t) = 



(3.50) 









26 2 ||B(T^7|jlD(u r ) - 2S 2 \\D(u r )\\]D(u r ) 




2 5 2 ||D(iC)jjE)(u r ) - 2<S 2 ||B(u r )|| B(u r ) 




25 2 \\Y>(%.)\\ B(u r ) - 2(5 2 ||B(u r )|| B(u r ) 



Since the stress tensors can be computed directly from the resolved field, (3.50) yields a 
time- and space-dependent formula for Cs(x, t). 

Thus, the DS-POD-ROM increases the viscosity coefficient by 

v DS :=2(C s (x 1 t)5) 2 \\D(u r )\l (3.51) 

where Cs(x, t) is the coefficient in ( |3.50 1, S is the lengthscale defined in § 3.2 and ||B(u r )|| 
the Frobenius norm of the deformation tensor B(u r ). Thus, the Dynamic Subgrid-Scale 
POD reduced-order model (DS-POD-ROM) has the form ( |3.11[ ), where 

6* (a) = -2<5 2 fv^.CKx,*) ||B(u r )|| VU+ o VUr N ) , (3.52) 



A km (a) = -2S 2 V^,Q(x,i)||D(u r 



(3.53) 



Note that v D s defined in (3.51) can take negative values. This can be interpreted as 
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backscatter, the inverse transfer of energy from high index POD modes to low index ones. 
The notion of backscatter, well-established in LES (see for example Sagaut 2006), was 
also found in a POD setting in the numerical investigation in Couplet et al. (2003). 



4. Numerical tests 

In this section, we use a structurally dominated 3D turbulent flow problem to test 
the four PO D-ROMs de scribed in §[3) (i) the ML-POD- ROM |3~14 | - pT5| ; (ii) the S- 
POD-ROM pl7|>-(|08l>; (ii i) the new VMS-POD-ROM p^ - pl^ ; and (iv) the new 
DS-POD-ROM pjgntjgl). We also include results for the POD-G-ROM K]\ (i.e., a 



POD-ROM without any closure model). A successful POD closure model should at least 



perform better than the POD-G-ROM (2.7). Finally, a DNS projection of the evolution 
of the POD modes served as benchmark for our numerical simulations: The closeness to 
the DNS data was used as a criterion for the success of the POD closure model. The 
qualitative behavior of all POD-ROMs is judged according to the following two criteria: 
(i) the kinetic energy spectrum, which represents the temporal average behavior of the 
POD-ROMs; and (ii) the time evolution of the POD coefficients, which measures the 
instantaneous behavior of the POD-ROMs. In § |4.1[ details of the numerical methods 
and parameter choices are given. In §4.2 numerical results are presented and discussed. 



4.1. Numerical Methods and Parameter Choices 

We investigate all four POD-ROMs in the numerical simulation of 3D flow past a circular 
cylinder at Re = 1, 000. The wake of the flow is fully turbulent. A parallel fluid flow solver 
is employed on a 144 x 192 x 16 finite volume mesh on the time interval [0, 300] to generate 
the DNS data. For details on numerical discretization, the reader is referred to Appendix 
A in |Wang et oT1 ( |2011| ). 

Collecting 1, 000 snapshots of the velocity field (m, U2, U3) over the time interval [0, 75] 



and applying the method of snapshots developed in |Sirovich ( 1987), we obtain the POD 
basis {<fi, • • • , <Pn}- These POD modes are then interpolated onto a structured quadratic 
finite clement mesh with nodes coinciding with the nodes used in the original DNS finite 
volume discretization. The first r = 6 POD modes capture 84% the system's kinetic 
energy. These modes are used in all POD-ROMs that we investigate next. For all the 
POD-ROMs, the time discretization was effected by using the explicit Euler method with 
At = 7.5 x 10~ 4 . 



It is important to note that the quadratic nonlinearity in the NSE (2.1 ) allows for easy 



precomputation of the vector b, the matrix A and the tensor B in the POD-G-ROM ( 2.8 ) 



For the general nonlinear EV POD closure model (3.11), however, the vector b(a) and 



the matrix A (a) that correspond to the additional closure terms have to be recomputed 
(reassembled) at each time step. Since the POD basis functions are global, although only 
a few are used in POD-ROMs (r <C N), reassembling b(a) and A(a) at each time step 
would dramatically increase the CPU time of the corresponding POD-ROM. Thus, the 
major advantage of POD-ROMs (the dramatic decrease of computational time), would 
be completely lost. 

To ensure a high computational efficiency of the POD-ROMs, we utilize two ap- 
proaches: (1) Instead of updating the closure terms in the POD-ROMs every time step, 
we recompute them every 1.5 time units (i.e., every 20,000 time steps). The previous 
numerical investigations in Wang et al. ( |2011 ) showed that this approach does not com- 
promise the numerical accuracy of the S-POD-ROM ( 3.17 )-( 3.18 ). (2) We employ the 



two-level algorithm introduced in Wang et al. (2011) to discretize the nonlinear closure 
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models. Before briefly describing the two-level algorithm, we emphasize that, in order to 
maintain a fair numerical comparison of the four POD-ROMs, we used both algorith- 
mic choices listed above in all four POD-ROMs. Therefore, the success or failure of the 
POD-ROM can solely be attributed to the closure term, which is the only distinguishing 
feature among all POD-ROMs, and not to the specific algorithmic choices, which are the 
same for all POD-ROMs. 

The two-level algorithm used in all four POD-ROMs is summarized below. 



£ = 0; compute b, A,B on the fine mesh ; 
for I = to M - 1 

compute b(a ), A(a ) on the coarse mesh 

a i+i := j?(a*); 
endf or 



(4.1) 



two-level 
algorithm 



In (4.1), M represents the total number of time steps. The idea in the two-level algo- 



rithm is straightforward: Instead of computing the closure terms b(a ), A(a ) directly 
on the fine mesh (as done in the standard one-level algorithm), the two- level algorithm 
discretizes them on a coarser mesh. Thus, the two-level algorithm is much more effi- 



cient than the standard one- level algorithm. Indeed, in Wang et al. (20111 it was shown 



that the two-level algorithm (4.1) achieves the same level of accuracy as the one-level 



algorithm while decreasing the computational cost by an order of magnitude. In all four 
POD-ROMs, we apply the two-level algorithm with a coarsening factor R c = 4 in both 
radial and azimuthal directions. Thus, the vectors and matrices related to the nonlinear 
closure terms are computed on a coarse finite element mesh with 37 x 49 x 17 grid points. 
In §|3.2| we proposed two definitions for the POD lengthscale S. Since in the finite 



element discretization that we employ, definition (3.10) is harder to implement than 



(3.9), we use the latter. Thus, using definition (|3.9|) with r — 6, we obtain 6 — 0.1179, 



which is the POD lengthscale that we will use in all four POD-ROMs. For the DS-POD- 



ROM (3.52 )-(3.53 1, we need to define the second (test) Galerkin projection (3.41) and 
the corresponding filter radius S. Choosing R = 1 in (3.41) and using (3.9), we obtain 
5 = 0.1769. 

The constants in EV LES models are determined in a straightforward fashion, utiliz- 



ing scaling laws satisfied by general 3D turbulent flows (see for example Sagaut 2006) 



Although the energy cascade concept in a POD context was verified numerically in Cou-| 
plet et al. (2003), there are no general scaling laws available in this setting. Thus, to our 



knowledge, the correct values for the EV constants a in the ML-POD-ROM ( 3. 14 )-( 3.15 1 



andCs in the S-POD-ROM ( 3.17 )-( 3.18) and the new VMS-POD-ROM ( 3.25 )-( 3.34) are 



still not known. To determine these EV constants, we run the corresponding POD-ROM 
on the short time interval [0, 15] with several different values for the EV constants and 
choose the value that yields the results that are closest to the DNS results. This approach 
yields the following values for the EV constants: a = 3 x 10~ 3 for the ML-POD-ROM, 
C s = 0.1426 for the S-POD-ROM, and C s = 0.1897 for the VMS-POD-ROM. We em- 
phasize that these EV constant values are optimal only on the short time interval tested, 
and they might actually be non-optimal on the entire time interval [0, 300] on which 
the POD-ROMs are tested. Thus, this heuristic procedure ensures some fairness in the 
numerical comparison of the four POD-ROMs. 

In the VMS-POD-ROM, only the first POD basis is considered as the large resolved 
POD mode, that is, rj, = 1 in (3.20 ). In the DS-POD-ROM, since vns can be negative, we 



use a standard "clipping" procedure to ensure the numerical stability of the discretization 
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(see for example Sagaut||2006 ) . Specifically, we let Cs(x,t) = max{Cs(x, t), — 0.2}. The 



value —0.2 is determined as follows: We first run the DS-POD-ROM without "clipping" 
for the time interval [0, 15] and record Cg , the average negative value of Cs(x,t). We 
then run on the entire time interval [0, 300] the DS-POD-ROM with a "clipping" value 
Cg ave /2 = —0.2. We note that there are alternative procedures to deal with the same 



issue in LES, such as VDSMwc (Morinishi & Vasilyev 2002 1. We utilized the standard 



"clipping" procedure described above as a first step in the numerical investigation of the 
DS-POD-ROM. 

4.2. Numerical results 
In th is secti on, we test the four P OD-R O Ms d escribed in §[3j (i) the ML-POD -ROM 



(3.14)-(3. 151; (ii) the S-POD-ROM ( 3.17 )-( 3.18|; (iii) the new VMS-POD-ROM (3.25) 



(3.34); and (iv) the new DS-POD-ROM (3.52)-(3.53). To assess their performance, we 



compare these four POD-ROMs with the POD-G-ROM §2J] (i.e., POD-ROM without 
any closure model) and the DNS projection of the evolution of the POD modes. A 
successful POD-ROM should perform significantly better than the POD-G-ROM and 
yield results that are close to those from the DNS. The POD-ROM numerical results 
are judged according to the following two criteria: (i) the kinetic energy spectrum, which 
represents the temporal average behavior of the POD-ROMs; and (ii) the time evolution 
of the POD coefficients, which measures the instantaneous behavior of the POD-ROMs. 
We also include a computational efficiency assessment for all four POD-ROMs. 

Before starting the quantitative comparison of the four POD-ROMs, we first give a 
flavor of the topology of the resulting flow fields. Figure[T]presents snapshots of horizontal 
velocity at t = 142.4 s for DNS, POD-G-ROM, ML-POD-ROM, S-POD-ROM, VMS- 
POD-ROM, and DS-POD-ROM. For clarity, only five isosurfaces are drawn. Taking the 
DNS results as a benchmark, the POD-G-ROM seems to add unphysical structures. 
The ML-POD-ROM, on the other hand, appears to add too much numerical dissipation 
to the system and thus eliminates some of the vortical structures in the wake. The 
S-POD-ROM, VMS-POD-ROM, and DS-POD-ROM perform well, capturing a similar 
amount of structure as the DNS. It also seems that there is some phase shift for all these 
POD-ROMs. Due to space limitations, only one time instance snapshot is shown for the 
POD-ROMs. The general behavior over the entire time interval is similar; it can be found 
at http : / /www. math . vt . edu/people/wangzhu/P0D_3DNumComp .html. 

Figure [2] presents the energy spectra of the four POD-ROMs and, for comparison 
purposes, of the POD-G-ROM. The five energy spectra are compared with the DNS 
energy spectrum. All energy spectra are calculated from the average kinetic energy of the 
nodes in the cube with side 0.1 centered at the probe (0.9992, 0.3575, 1.0625). It is clear 
that the energy spectrum of the POD-G-ROM overestimates the energy spectrum of the 
DNS. The energy spectrum of the ML-POD-ROM, on the other hand, underestimates the 
the energy spectrum of the DNS, especially at the higher frequencies. The S-POD-ROM 
has a more accurate spectrum than the ML-POD-ROM, but it displays high oscillations 
at the higher frequencies. The VMS-POD-ROM is a clear improvement over the S-POD- 
ROM, with smaller oscillations at the higher frequencies. The energy spectrum of the DS- 
POD-ROM is qualitatively similar to that of the VMS-POD-ROM. The DS-POD-ROM 
spectrum decreases the amplitude of the high frequency oscillations of the VMS-POD- 
ROM even further, although it introduces some sporadic large amplitude oscillations 
at high frequencies. To summarize, the DS-POD-ROM and the VMS-POD-ROM yield 
the most accurate energy spectra, i.e., the closest to the DNS energy spectrum. On the 
average, the DS-POD-ROM performs slightly better than the VMS-POD-ROM. 

As the second criterion in judging the performance of the four POD-ROMs, the time 
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Figure 1. (Continued on n ext p age.) Snapshots of horizo ntal v e locity at t — 142.5 s for: (a) 



3.17 




3.18 









(e) the new VMS-POD-ROM p35|-([334||; and (Tpthe new DS-POD-ROM 



Five isosurfaces are plotted. Note that the POD-G-ROM adds unphysical struc- 



VMS-POD-ROM, and DS-POD-ROM perform well, capturing a similar amount of structure as 
the DNS. 



(a) 




DNS 
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(f) 




DS-POD-ROM 



evolutions of the POD basis coefficients ai(-) and a^-) on the entire time interval [0, 300] 
are shown in Figures 3]|4 We note that the other POD coefficients have a similar behavior. 
Thus, for clarity of exposition, we include only ai(-) and 04(-). The POD-G-ROM's time 
evolutions of a\ and 0,4 are clearly inaccurate. Indeed, the magnitude of 04 is nine times 
larger than that of the DNS projection, which indicates the need for closure modeling. The 
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ML-POD-ROM's time evolutions of a\ and 04 are also inaccurate. Specifically, although 
the time evolution at the beginning of the simulation (where the EV constant a was 
chosen) is relatively accurate, the accuracy significantly degrades toward the end of the 
simulation. For example, the magnitude of 04 at the end of the simulation is only one 
eighth of that of the DNS. The S-POD-ROM yields more accurate time evolutions than 
the ML-POD-ROM for both a\ and 04, although the magnitude of the POD coefficients 
stays almost constant at a high level. The VMS-POD-ROM's time evolutions of ai and 04 
are better than those of the S-POD-ROM, since the magnitudes of the POD coefficients 
are closer to those of the DNS. Finally, the DS-POD-ROM also yields accurate results. We 
note that the DS-POD-ROM's a\ and 04 coefficients have significantly more variability 
than the corresponding coefficients of the VMS-POD-ROM. This is a consequence of 
the fact that the EV coefficient Cs varies in time and space for the DS-POD-ROM, 
whereas it is constant for the VMS-POD-ROM. To summarize, the DS-POD-ROM and 
the VMS-POD-ROM perform the best. On the average, the DS-POD-ROM performs 
slightly better than the VMS-POD-ROM. 

Based on the energy spectra and the the time evolutions of the POD basis coefficients 
ai(-) and a 4 (-), the DS-POD-ROM and the VMS-POD-ROM consistently outperform the 
ML-POD-ROM and the S-POD-ROM. To determine which one of the DS-POD-ROM 
and the VMS-POD-ROM performs best, we collected the results in Figures |4]Jd)-[4]Je) 
(corresponding to the time evolution of the POD basis coefficient (n(-) for the DNS 
projection, the VMS-POD-ROM and the DS-POD-ROM) and we displayed them in the 
same plot in Figure [5] Since it is difficult to distinguish between the results from the 
VMS-POD-ROM and the DS-POD-ROM, we zoomed in on the POD basis coefficient 
04 over the time interval [266,282]. Based on the plot in the inset, it is clear that, for 
this time interval, the DS-POD-ROM performs better than the VMS-POD-ROM. More 
importantly, it appears that the magnitude of a± in the DS-POD-ROM displays some 
of the variability displayed by the DNS; the magnitude of the VMS-POD-ROM's a 4 
coefficient, on the other hand, displays an almost periodic behavior. We believe that the 
variation of the DS-POD-ROM's 04 coefficient is due to the dynamical computation of 
the EV coefficient, which changes both in space and time; the EV coefficient of the VMS- 
POD-ROM, however, is constant and is computed at the beginning of the simulation. 

To gain further insight into the behavior of the DS-POD-ROM and the VMS-POD- 
ROM, we considered the root mean square horizontal velocity of the two models. Figure[6] 
presents the average horizontal velocity (u) , which is computed at points with coordinates 

x = 3.2937 and y = 2.2796, and the root mean square horizontal velocity u rms = 
1/2 

(u — (u) , u — (u)) 1 I (u), where (•) here denotes the spatial average in the z-direction. 
Both the DS-POD-ROM and the VMS-POD-ROM yield average horizontal and root 
mean square velocities that are in close agreement with the DNS data. As for the time 
evolution of the POD basis coefficient a 4 in Figure [5j the DS-POD-ROM and the VMS- 
POD-ROM results are practically indistinguishable. 

To summarize, the VMS and DS approaches, which are state-of-the-art closure model- 
ing strategies in LES, yield the most accurate POD closure models for the 3D turbulent 
flow that we investigated. A natural question, however, is whether the new POD closure 
modeling strategies that we proposed display a high level of computational efficiency, 
which is one of the trademarks of a successful POD-ROM. To answer this question, we 
computed the CPU times for all four POD- ROMs and compared them with those of the 
DNS and the POD-G-ROM. 

To make such a comparison, however, we first need to address the numerical differ- 
ences between the DNS and the POD-ROMs. First, the discretizations used in the two 
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Table 1. Speed-up factors of POD-ROMs. 



POD-G-ROM 


ML-POD-ROM 


S-POD-ROM 


VMS-POD-ROM 


DS-POD-ROM 


S f 665 


659 


36 


41 


23 



approaches are completely different. Indeed, the spatial discretization used in the DNS 
was the finite volume method, whereas for the POD-ROMs we used a finite element 
method. Furthermore, the time-discretization used in the DNS was second-order (Crank- 
Nicolson and Adams-Bashforth) , whereas in the POD-ROMs we used a first-order time 
discretization (explicit Euler). The time steps employed were also different: At = 2x 1CU 3 
in the DNS and At = 7.5 x 1(T 4 in the POD-ROM. Most importantly, the DNS was 
performed on a parallel machine (on 16 processors), whereas all the POD-ROM runs were 
carried out on a single-processor machine. Thus, to ensure a more realistic comparison 
between the CPU times of the DNS and the POD-ROMs, we multiplied the CPU time 
of the DNS by a factor of 16. 

To measure the computational efficiency of the four POD-ROMs, we define the speed- 
up factor 

CPU time of DNS 

f ~ CPU time of POD-ROM ^ ' ' 

and list results in Table [T] The most efficient model is the POD-G-ROM. This is not 
surprising, since no closure model is used in POD-G-ROM and thus no CPU time is spent 
computing an additional nonlinear term at each time step. The second most efficient 
model is the ML-POD-ROM. This is again natural, since only a linear closure model is 
employed in the ML-POD-ROM and thus the computational overhead is minimal. The 
speed-up factors for the S-POD-ROM, the VMS-POD-ROM and the DS-POD-ROM are 
one order of magnitude lower than those for the ML-POD-ROM and the POD-G-ROM. 
The reason is that the former use nonlinear closure models, which increase significantly 
the computational time. Note, however, that the S-POD-ROM, the VMS-POD-ROM and 
the DS-POD-ROM are still significantly more efficient than the DNS. 
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Figure 2. Ki neti c energy spectrum of the DNS (blue) and the POD-RO Ms (re d ): (a) the 
POD-G-ROM jO overestimates the DNS spectrum; (b) the ML- POD-ROM | |3.14| |- |3.15| un- 



derestimates the DNS spectrum; (c) the S-POD-ROM ( |3.17 1- ( 3. 18| ) yields a more accurate spec- 
trum than the ML-PO D-ROM, but displays high oscillations at the higher frequencies; (d) the 



new VMS-POD-ROM ( |3.25| )-( |3.34| clearly improves the accuracy of the S-POD-ROM's spec- 
trum, displa ying smaller oscillations at the higher frequencies; and (e) the new DS-POD-ROM 
( 3.52 H 3.53 1 decreases the amplitude of the high oscillations of the VMS-POD-ROM's spectrum, 
although it displays sporadic undershoots. 
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Figure 3. Time evolutions of the POD b asis coefficient a\ of the DNS (blue) and 
the POD-ROMs (red ): (a) the POD-G-ROM pjj overestimates the DNS res ults; ( b ) the 
ML-POD-ROM pU4} -< pU5l l underestimates the DNS results; (c) the S-POD-ROM (f3Tl7b - d37T8 > 
yields more accurate resultsthan the ML-POD-ROM; (d) the new VMS-POD-ROM |3"35) -( pHI l 
improves the accuracy of the S-PO D-ROM res ults, especially toward the end of the simulation; 
and (e) the new DS-POD-ROM ( 3.52 1- ( 3.53 1 yields results that are similar to those of the 
VMS-POD-ROM. 
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Figure 4. Time evolutions of the POD b asis coefficient 124 of the DNS (blue) and 



ML-POD-ROM 1 3.14 1-13.15 1 underestimates the DNS results; (c) the S-POD-ROM 



yields more accurate resuTtsthan the ML-POD-ROM; (d) the new VMS-POD-ROM _ 
improves the accuracy of the S- POD-ROM results, especially toward the end of the simulation; 



3.17 




3.18 






3"3l 



and (e) the new DS-POD-ROM | |3.52| -| |3.53 1 performs slightly better than the VMS-POD-ROM, 
recovering some of the variability of the DNS results. 
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5. Conclusions 

This paper put forth two new POD-ROMs (the DS-POD-ROM and the VMS-POD- 
ROM), which are inspired from state-of-the-art LES closure modeling strategies. These 
two new POD-ROMs together with the ML-POD-ROM and the S-POD-ROM were tested 
in the numerical simulation of a 3D turbulent flow past a cylinder at Re = 1, 000. For 
completeness, we also included results with the POD-G-ROM (i.e., a POD-ROM without 
any closure model), as well as the DNS projection of the evolution of the POD modes, 
which served as benchmark for our numerical simulations. 

To assess the performance of the POD-ROMs, two criteria were considered in this 
paper: the kinetic energy spectrum and the time evolution of the POD basis coefficients. 
The former is used to measure the average behavior of the POD-ROMs and the latter is 
used to quantify the instantaneous behavior of these models. Both the POD-G-ROM and 
the ML-POD-ROM yielded inaccurate results. The DS-POD-ROM and the VMS-POD- 
ROM clearly outperformed these two models, yielding more accurate results for both 
the kinetic energy spectrum and the time evolution of the POD basis coefficients. The 
DS-POD-ROM performed slightly better than the VMS-POD-ROM for both criteria and 
also seemed to display more adaptivity in terms of adjusting the magnitude of the POD 
basis coefficients. Overall, however, the two models yielded similar qualitative results. 
This seems to reflect the LES setting, where both the DS and the VMS closure modeling 
strategies are considered state-of-the-art (Hughes et al. 2001a|fe ). The DS-POD-ROM and 



the VMS-POD-ROM, although not as computationally efficient as the POD-G-ROM or 
the ML-POD-ROM, significantly decreased the CPU time of the DNS. To summarize, 
for the 3D turbulent flow that we investigated, the DS-POD-ROM and the VMS-POD- 
ROM were found to perform the best among the POD-ROMs investigated, combining a 
relative high numerical accuracy with a high level of computational efficiency. 

We plan to further investigate several other research avenues. First, we plan to study 
more efficient time-discretization approaches and take advantage of parallel computing 
in order to further decrease the computational time and, at the same time, increase the 
dimension (and the thus physical accuracy) of the POD-ROMs. Second, since the linear 
closure model (ML-POD-ROM) is computationally efficient, but only works on a relative 
short time interval if the appropriate EV coefficient a is chosen, we will investigate a 
hybrid approach: We will use the DS approach to calculate a only when the flow displays 
a high level of variability, and then use this value in the ML-POD-ROM as long as the flow 
does not experience sudden transitions. Third, using these computational developments, 
we will investigate the new POD-ROMs in more challenging, higher Reynolds number 
structurally dominated turbulent flows. Finally, we plan to employ the new POD-ROMs 
in other scientific and engineering applications in which accurate POD closure modeling 
is needed, such as optimal control, optimization, and data assimilation problems. 

We greatly appreciate the financial support of the Air Force Office of Scientific Research 
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